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We present results from Monte Carlo calculations investigating the properties of the homogeneous, 
spin-balanced unitary Fermi gas in three dimensions. The temperature is varied across the superfiuid 
transition allowing us to determine the temperature dependence of the chemical potential, the energy 
per particle and the contact density. Numerical artifacts due to finite volume and discretization are 
systematically studied, estimated, and reduced. 
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I. INTRODUCTION 

Ultracold, dilute gases of fermionic atoms have been 
studied extensively lately, in part due to the system be¬ 
ing the simplest environment with strong interactions be¬ 
tween fermions (for recent reviews of this active field see 
e.g. Refs. [Ilia). Most remarkably, a three-dimensional 
(3d) atomic gas with two hyperfine states, say lithium or 
potassium, can be constructed to have resonant interac¬ 
tions: by applying an external magnetic field the S-wave 
scattering length a can be tuned to satisfy 1/a = 0. In 
this special situation, the properties of the gas become 
universal, dependent only on the density and tempera¬ 
ture. What can be learned in the atomic physics labora¬ 
tory then has implications for nonrelativistic fermions as 
small as nucleons. 

This resonant Fermi gas - often called the “unitary 
Fermi gas” due to the scattering being limited only by 
unitarity - provides an excellent opportunity for quan¬ 
titative theoretical calculations. The clean separation of 
scales means that much can be inferred from dimensional 
analysis and scaling arguments. What remains to be de¬ 
termined are universal dimensionless constants, such as 
the Bertsch parameter [5], or functions of the product 
of the inverse temperature /? and the chemical potential 
M i which completely specify the thermodynamic and 
hydrodynamic behavior of the unitary Fermi gas. Much 
effort has gone into using first-principles numerical meth¬ 
ods to determine these quantities (see Refs. m\ for re¬ 
views) . 

In this paper we use lattice Monte Carlo methods to 
give numerical results for thermodynamic quantities as 
the temperature is varied through the superfluid phase 
transition. In particular we determine the chemical po¬ 
tential, mean energy density, and contact density as func¬ 
tions of temperature. The corresponding universal func¬ 
tions /(/3^) are made dimensionless by taking ratios with 
the appropriate powers of the Fermi energy sp. It is no¬ 
table that several other numerical methods for studying 
the Fermi gas cannot study the superfluid phase. With 
this study, we also pay particular attention to investi¬ 


gating and quantifying the systematic uncertainties as¬ 
sociated with taking the thermodynamic limit and the 
continuum limit, which is crucial to obtain correct physi¬ 
cal results Generally, we find good agreement with 

experiment. 

In previous work m we used the Diagrammatic De¬ 
terminant Monte Carlo (DDMC) algorithm [TTHI^ to 
numerically determine the critical temperature Tc and 
thermodynamic properties of the unitary Fermi gas at 
T = Tc- Here we study the temperature dependence 
of physical observables in the approximate range Tcj2 < 
T < 2Tc. An approach which is formulated in the contin¬ 
uum, bold-line diagrammatic Monte Carlo (bold DMC), 
has been used to compute quantities above the critical 
temperature [Hill], and these results agree well with 
experimental measurements. This method as presently 
formulated does not extend to temperatures below Tc 
due to the singularity in the pair propagator appearing in 
the superfluid phase. Temperature effects have also been 
studied in a lattice computation using a hybrid Monte 
Carlo approach |16j : however results using different lat¬ 
tice spacings are not presented there. 

II. SETUP 

We consider a system of equal-mass fermions with two 
spin components labeled by the spin index a = {t,)-}. 
Since the details of the physical potential governing the 
interatomic interactions are irrelevant in the dilute limit 
realized in cold-atom experiments, we can work on a spa¬ 
tial lattice provided that we also take the dilute limit [m. 
The Hamiltonian is that of the simple Fermi-Hubbard 
model in the grand canonical ensemble, 

H = ^^(ck ~ Ta)c\^cr^^^cr + U '"xtCx^, (1) 

k,cr X 

where the first term corresponds to the kinetic part of 
the Hamiltonian Hkin and the second term to the in¬ 
teraction part Flint- The units are chosen such that 
h = kp = 2m = 1. We work on a 3d periodic lattice 
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with lattice spacing b and sites. The discrete dis¬ 
persion relation reads Ck = — cos kjb); 

is the chemical potential and the fermionic creation 
operator. The coupling constant U < 0 corresponding to 
attractive interaction can be tuned so that the scattering 
length becomes infinite. The corresponding value in the 
infinite-volume limit is U = —7.914/6^ which is the value 
we use throughout the calculation. Another approach is 
to include finite-volume effects in this two-body match¬ 
ing calculation [MUSI ESI- It remains to be seen which 
approach leads to a milder extrapolation of many-body 
results to the continuum limit. 

The partition function Z = Trexp(—/377) can be writ¬ 
ten as a series of products of two matrix determinants 
built of free finite-temperature Green’s functions m- If 
~ k-l = Mi 9-8 is always assumed to be the case in the 
present work, the two determinants are identical since 
the spin-dependence enters only via the chemical poten¬ 
tial. Consequently all terms in the series are positive, 
and the series can be used as a probability distribution 
for Monte Carlo sampling. 

We use the DDMC algorithm as introduced in m with 
several modifications which increase the efficiency by re¬ 
ducing autocorrelation effects as compared to the original 
setup. We account for remaining autocorrelations by bin¬ 
ning the data to the point where the statistical error is 
insensitive to the bin size. A detailed description of our 
numerical setup is given in [MEniEI]. 

We performed many calculations, varying the dimen¬ 
sionless inputs for the chemical potential {^J.b'^) and in¬ 
verse temperature (/3/6^), as well as the number of lattice 
points L^, so that controlled extrapolations to the ther¬ 
modynamic and continuum limits could be taken for a 
range of temperatures in both the normal and superfluid 
phases. Each diamond in Fig. [^represents the thermo¬ 
dynamic limit of lattice calculations done for values of 
and /3/6^ resulting in the corresponding filling factor 
ly and dimensionless product /3/i. 


III. THERMODYNAMIC AND CONTINUUM 
LIMITS 

We set the physical scale via ly = nb^, where ly = 
'^xo'Cxcr) is the dimensionless hlling factor and n the 
particle number density. Due to universality all physical 
quantities are given in units which can be expressed as 
appropriate powers of the Fermi energy ep = 

To extract the physical results we need to perform two 
limits: the thermodynamic limit to infinite system size 
and then the continuum limit to zero lattice spacing. 

First we take the thermodynamic limit in order to es¬ 
timate and reduce systematic errors due to finite volume. 
For each {iib'^,/3/b'^) we perform computations with sev¬ 
eral (usually 3 or 4) sizes of cubic volumes, V = L^, 
and extrapolate results as IjL ^ 0. At the lowest fill¬ 
ing factor (smallest b) we used volumes up to U = 32^ 
(corresponding to up to 1000 particles). Typical ranges 
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FIG. 1. Diamonds correspond to parameters used for indi¬ 
vidual simulations extrapolated to the thermodynamic limit. 
The 15 bands indicate the subsets of data included in the 
continuum (ty —>■ 0) extrapolations as described in Sec. |III[ 


are L = 6 to L = 14 at higher filling factors (about 40 
to 700 particles) and L = 10 to L = 32 at lower filling 
factors (about 100 to 1000 particles). We find that the 
filling factor ly is the quantity most sensitive to finite- 
volume effects. With periodic boundary conditions, we 
expect finite-volume effects to cause an increase in filling 
factor compared to the infinite-volume limit. In small 
volumes the particles will feel an enhanced attraction 
not just from their neighboring particles, but also from 
their round-the-world doppelgangers. In agreement with 
m we observe that the data for the hlling factor are 
Ht well by a linear function of 1/L. Once we have ly 
in the thermodynamic limit, we can obtain ep and the 
dimensionless observables by taking v to the appropri¬ 
ate power and multiplying by quantities which have no 
statistically signihcant hnite-volume errors. Within the 
statistical uncertainties, it appears that hnite-volume er¬ 
rors are negligible for E/N and Eint/N. Examples of 
thermodynamic extrapolations are shown in Figs. and 
[^ for different parameter sets in the superhuid and nor- 
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PH=2.7, p/b^=9, nb^=0.3 
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FIG. 2. Examples of thermodynamic limit extrapolations in the snperfluid phase for the filling factor v (left), the energy per 
particle EjN = EjjJ'v (middle), and the interaction energy per particle Ei-ntlN = Eint/L^v (right), which yields the valne of 
the contact. 


mal phases, respectively. 

For the continuum limit we vary the dimensionless 
chemical potential fib"^ such that the filling factor tends 
to zero. This is equivalent to 6 —)■ 0 since b oc if 
n is fixed to be a constant, physical value. Dimension¬ 
less ratios of physical quantities can then be extrapolated 
to the continuum limit by assuming discretization effects 
can be parametrized using a power series in b, or equiv¬ 
alently 

Our previous work |10[ 122) focused on determining the 
critical temperature and computing thermodynamic ob¬ 
servables there. The critical temperature was determined 
at several values of the lattice spacing. In practice this 
was done by varying the inverse temperature and chem¬ 
ical potential in lattice units to the point where a finite- 
size scaling study of the order parameter indicated the 
transition would occur in infinite volume. Here we extend 
that work by taking the continuum limit of observables 
for a range of temperatures on either side of the phase 
transition. 


Ideally the continuum limit would be taken varying b 
along lines of constant The numerical data acquired, 
however, do not lie exactly along lines of constant /3^. 
Therefore, we group the data in several narrow bands of 
y = /3/i values and extrapolate the data within each band 
to the continuum limit. Each band is defined by a cen¬ 
tral value ?/o Emd a width, as shown in Fig. In order 
to account for mild y dependence, we find it sufficient 
to introduce a term proportional to Sy = y — yo in the 
extrapolation function. The lattice spacing dependence 
of a physical quantity X can be written as a power se¬ 
ries in the lattice spacing b oc [53]. Therefore, our 
continuum limit fits are to functions of the form 

X{yo-,v,5y)= Xo{^ + di5y + ^Ckv''^^ ■ (2) 

The fit parameter Xq = Xo(j/o) is then taken to be the 
continuum limit result. The other fit parameters, di and 
Cfc also depend on yo, but we suppress this dependence 
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PH=2.21, p/b^=2.764, nb^=0.8 


PH=2.21, p/b®=2.764, nb^=0.8 
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FIG. 3. Examples of thermodynamic limit extrapolations in the normal phase for the filling factor v (left), the energy per 
particle EjN — E/L^v (middle), and the interaction energy per particle Eint/N = EintfL^i' (right), which yields the value of 
the contact. 


in the notation. In almost every case the Monte Carlo 
data are sufficient to determine ci but not the coeffi¬ 
cients of higher-order terms. In other words, the data 
points indeed look linear in However, given that, 

especially for larger y, the numerical values of are 
not very small, it is prudent to allow for higher-order con¬ 
tributions in the numerical data. Therefore we introduce 
Bayesian prior distributions for the Ck with fc > 1 |24j . 
In the cases where the Monte Carlo results do constrain 
C 2 (i.e. fits to yje-p at low y) we find C 2 « 0.3. There¬ 
fore, we take Gaussian prior distributions centered at 0 
with width 0.3. We found very little difference in the 
fits where we set K = 2 or K — ?>, but we used the lat¬ 
ter, more conservative, option for the results presented 
here. Finally, we also performed fits which included a 
term fi 6y with a Gaussian prior for fi of 0.0 ± 1.0. 
This had no effect on the fits, so for simplicity we omit 
this term from our final fits. 

In Fig. 1^ we show the results of these fits for 5 of the 15 
bands. The fit curve and corresponding error are evalu¬ 


ated at the value of /3/r given in the legend, and represent 
an interpolation of the data to the central value of the 
lettered band shown in Fig. as well as the continuum 
extrapolation to z/ = 0. In the case of band F, we have 
several data points generated with /3/i = 2; we emphasize 
these points with stars. Within uncertainties, the widths 
of the bands in jSy are sufficiently narrow that interpo¬ 
lating in Py is mild, well-parametrized by the di term in 
Eq. 

Summarizing the size of discretization error for the 
continuum extrapolations, we show the values obtained 
for the coefficient Ci of in Eq. ([^. Figure shows 
that the Monte Carlo data for the chemical potential and 
energy per particle have significant linear dependence on 
the lattice spacing, as parametrized by especially 
at lower Py. This is similar to what was seen for lat¬ 
tice determinations of Tc [T0l[T2]. Nevertheless, the fits 
described above, which include possible contributions of 
higher-order terms, include estimates of these lattice ar¬ 
tifacts in the uncertainties quoted below. We tried a 
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FIG. 4. Continuum limit extrapolation along bands of constant /3/r±(5(/3p) (see Fig.[l| for the chemical potential (left), energy 
density (middle) and contact density (right). See further discussion at the end of Sec. |III| regarding the interpolation in /3/r. 
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stood by looking at the lattice Monte Carlo estimator for 
E\^i-a/N = [in], which can be expressed as 

-E'kin __ g A So'('^X(TC(x+j)cr) 

L^V I V 

This expression contains the ratio of the kinetic energy 
operator X)cr(‘'xcrC(x+j)(T) filling factor operator 

^^(c^^Cxcr)- These two operators have a very simi¬ 
lar structure, which explains to some extent why finite- 
volume errors are negligible within the error bars of the 
data. The same holds for the interaction part of the en¬ 
ergy. Therefore it is sufficient to consider the finite-size 
scaling of l/ep (which follows directly from the finite-size 
scaling of z/), while the data for ElLi^v obtained at dif¬ 
ferent lattice sizes can be fitted by a constant. Our data 
confirms this scaling, as the constant fits of E'/L^iz yield 
acceptable y^-values, see the middle panel of Figs. and 
[^for several examples. 

The results for the energy per particle E/Epc, where 
= (3/5)iVeF is the ground state energy of the free 
gas, are shown in the right panel of Fig. Like for the 
chemical potential, we obtain excellent agreement with 
experimental data US] US] and theory [I1I27I- 



C. Contact density 


FIG. 5. The fit parameter Ci obtained from continuum ex¬ 
trapolations, see Eq. Q, of Monte Carlo data for the chemi¬ 
cal potential (top), energy per particle (middle), and contact 
density (bottom). 

variety of other fits, altering the bands in /?/r, and omit¬ 
ting data with large filling factor, e.g., > 0.5 [9|. 

These variations produced fits which were in agreement 
with our final results, but were less precise in some cases. 

IV. RESULTS 

A. Chemical potential 

The left panel of Fig. shows the continuum limit 
of the chemical potential as a function of /3/r. We see 
excellent agreement with experimental data [551 [55] , as 
well as with several other theoretical predictions [HUT]. 
Our results below Tc capture the experimentally observed 
change of the slope of the chemical potential curve. 

B. Energy per particle 

The energy is composed of the kinetic energy E]^in and 
the interaction energy Flint. We find that, within the 
statistical uncertainties, neither Ey^i^/N nor Ei^t/N ex¬ 
hibit dependence on L. This insensitivity can be under- 


The contact density can be interpreted as a measure 
of the local pair density |28| . The contact plays a crucial 
role for several universal relations derived by Tan Ei- 
I3II- We use the definition C = m^goEint, where go is 
the physical coupling constant |28l [32] . The contact is 
related to the contact density C via C = J C(r)d^r, or for 
homogeneous systems simply C = CV. 

In E2] we have presented results for the contact density 
at the critical point. Now we extend this study to other 
values of the temperature. For the finite-size scaling we 
can rewrite the dimensionless contact density as 

C _ UEi^y ^ U Flint V Flint /.N 

ef “ 4F34 (37r2)4/3 

Since Ei^t/N is independent of L within uncertainties, 
this part of the contact density for different lattice sizes 
is fit to a constant (see the right panel of Figs. andfor 
examples of such fits), while the thermodynamic limit for 
the part proportional to v~^l^ follows from the finite-size 
scaling of the filling factor v. 

Figure shows the contact density in the continuum 
limit. There has been recent progress experimentally in¬ 
vestigating Tan’s contact, mostly for trapped systems 
[551 1551 155] as well as numerical and analytical calcu¬ 
lations [m [51 iSHS]- The homogeneous contact in 
the normal phase has been studied experimentally in 
Ref. [33]. They find a sharp decrease in the contact 
around the superfluid phase transition. We do not ob¬ 
serve any such sudden change around Tc, but our results 
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FIG. 6. The chemical potential ^/ep (left panel) and the energy per particle E/Epq (right panel) in the continuum limit 
versus /3^. We compare our results (red circles; the empty circles denote our results at Tc from |10p with experimental data 
|25| (green squares), as well as results obtained with bold DMC [14] (blue triangles) and with the Luttinger-Ward formalism 
m (black dashed line). The gray bar indicates the critical point m with error margin. 
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FIG. 7. The contact density C/sp = C/kp in the continuum 
limit versus f5fi. We compare our results (red circles; the 
empty circle denotes our result at from |^) with experi¬ 
mental data [33] (green squares), as well as results obtained 
with bold DMC [15] (blue triangles), Luttinger-Ward formal¬ 
ism |34| (black dashed line) and the zero-temperature results 
from [35] (cyan line with error margin), [36] (purple dotted 
line with error margin) and |37| (orange dash-dotted line with 
error margin). 


above show good agreement with their data. Our re¬ 


sults at low temperature also show excellent agreement 
with the zero-temperature numerical m and experimen¬ 
tal results [SU |36] (the contact is not discussed explicitly 
in the latter reference, but can be easily extracted with 
the appropriate Tan relation yielding C/e| = 2C/57r = 
0.1184(64)). 

V. SUMMARY AND CONCLUSION 

In summary, we have calculated the chemical poten¬ 
tial, the energy density and the contact of a homogeneous 
3d balanced unitary Fermi gas at different temperatures 
above and below the critical point. Our results show good 
agreement with experimental measurements and provide 
a benchmark for future studies, in particular below Tc 
where few accurate predictions and measurements are 
available. 
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